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Abstract 

Vibrational energy relaxation (VER) of a selected mode in cytochrome c (hemeprotein) 
in vacuum is studied using two theoretical approaches: One is the equilibrium simulation 
approach with quantum correction factors, and the other is the reduced model approach 
which describes the protein as an ensemble of normal modes coupled with nonlinear coupling 
elements. Both methods result in estimates of VER time (sub ps) for a CD stretching 
mode in the protein at room temperature, that are in accord with the experimental data of 
Romesberg's group. The applicability of the two methods is examined through a discussion 
of the validity of Fermi's golden rule on which the two methods are based. 
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1 Introduction 

The harmonic (or normal mode) approximation has been a powerful tool for the analysis of few 
and many-body systems where the essential dynamics of the system consists of small oscillations 
about a well-defined mechanically stable structure. The concept of normal modes (NMs) is 
appealing in science because it provides a simple view for complex systems like solids and 
proteins. Though it had been believed that NMs may be too simplistic to analyze the dynamics 
of proteins, that is by no means always true: the experimental data of neutron scattering for 
proteins (B factor) indicate that the fluctuations for each residue are well represented by a 
simplified model using NMs It was also shown that such a large-amplitude motion as the 
hinge-bending motion in a protein is well described by a NM |2]. Importantly, NMs have been 
used to refine the x-ray structures of proteins Recently, large proteins or even protein 
complexes can be analyzed by using NMs El El • 

In this chapter, we are concerned with vibrational energy relaxation (VER) in a protein. 
This subject is related to our understanding of the functionality of proteins: At the most 
fundamental level, we must understand the energy flow (pathway) of an injected energy, that 
is channeled to do useful work. Due to the advance of larser technology, time-resolved spec- 
troscopy can detect such energy flow phenomena experimentally To interpret experimental 
data, and to suggest new experiments, theoretical approaches and simulations are essential 
as they can provide a detailed view of VER. However, VER in large molecules itself is still 
a challenging problem in molecular science 0. This is because VER is a typical many-body 
problem and estimations of quantum effects are difficult 9_ . There is a clear need to test and 
compare the validity of the existing theoretical methods. 

We here employ two different methods to estimate the VER rate in a protein, cytochrome c 
(see the next section for details). One is the classical equilibrium simulation method jlOj with 
quantum correction factors J2 El • The second is the reduced model approach ^Hl , which has 
been recently employed by Leitner's group jl4l I15j . The latter approach is based on NM con- 
cepts, which describes VER as energy transfers between NMs mediated by nonlinear resonance 
jl6j . We conclude with a discussion of the validity and applicability of such approaches. 



2 Cytochrome c 

Cytochrome c (cyt c) is one of the most thoroughly physicochemically characterized metallo- 
proteins ^1|^. It consists of a single polypeptide chain of 104 amino acid residues and is 
organized into a series of five a-helices and six /3-turns. The heme active site in cyt c consists 
of a 6-coordinate low-spin iron that binds Hisl8 and MetSO as the axial ligands. In addition, 
two cysteines (Cysl4 and Cysl7) are covalently bonded through thioether bridges to the heme 
(see Fig.^. Crystal structures of cyt c show that the heme group, which is located in a groove 
and almost completely buried inside the protein, is non-planar and somewhat distorted into 
a saddle-shape geometry. The reduced protein, ferrocytochrome c (ferrocyt c), is relatively 
compact and very stable, due to the fact that the heme group is neutral. 

The vibrational mode we have chosen for study is the isotopically labeled CD stretch in 
the terminal methyl group of the residue MetSO, which is covalently bonded to Fe in heme (see 
Fig. Our simulation model approximates the protein synthesized by Romesberg's group 
jl9j though their protein contains three deuteriums in MetSO (Met80-3D). The CH and CD 
stretching bands are located near 3000 cm~^ and 2200 cm~^, respectively. In contrast with 
the modeling of photolyzed CO in myoglobin jJUl, essentially a diatomic molecule in a protein 
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result, the modeling is more challenging: There is no clean separation between the system and 
bath modes because the CD bond is strongly connected to the environment. 




Figure 1: The structure of cytochrome c in the vicinity of the heme group, showing the 
thioether linkages and non-planar heme geometry. 



3 Quantum correction factor approach 

The classical Landau- Teller-Zwanzig theory of VER is attractive in that it allows us to base 
our estimate of the VER rate on a classical force autocorrelation function that contains the 
interaction coupling between the system and bath modes to all orders. The Hamiltonian for such 
a system is of the Caldeira-Leggett-Zwanzig form, where "bath" coordinates are represented 
as normal modes of the bath alone} The relaxing oscillator is introduced as a local "system" 
mode, coupled to the bath at all orders, including "bi-linear" coupling. 

Efforts have been made to introduce quantum effects through the use of "quantum correction 
factors." The dynamics of the classical system are computed, and the quantum effects are added 
a posteriori in a manner that accounts for the equilibrium quantum statistical distribution of 
the contributing quantum mechanical degrees of freedom. This approach is summarized below 
and applied to estimate the rate of VER for the CD bond in the terminal methyl group of Met 
80 in cytochrome c. 

3.1 Fermi's golden rule 

Our starting point for computing the rate of vibrational energy relaxation of the CD stretching 
mode in cytochrome c is Fermi's "golden rule" formula. The vibrational population relaxation 

^As shown below, the classical LTZ formula can be considered as a classical limit of the quantum mechanical 
population relaxation rate l/Ti. This result is derived by using both Fermi's golden rule and Bader-Berne 
theory ^D,. Though the transition rate itself can be derived without any assumption on the bath 

Hamiltonian, the Bader-Berne result stems from the assumption that the bath Hamiltonian is an ensemble of 
harmonic oscillators. 
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- dt COs{uJst) Cqm{t) = 77> Cqmi^s) (1) 



Ti phuJs/2 Jo ' ' pnujs/2 

where the force- force correlation function Cqm(i) is defined as 



Cqm(t) = ^iHt)m + mm)rira, (2) 



its Fourier transform is Cqm(i^), ^{t) is the quantum mechanical force applied to the system 
mode considered, ms is the system (reduced) mass, us is the system frequency, /3 is an inverse 
temperature, and the above bracket means a quantum mechanical average. Note that in the 
classical limit ^ — >• 0, the prefactor in front of the integral in Eq. becomes unity, and the 
expression reduces to the well-known classical VER formula. The issue is that this limit does 
not represent well the VER for high frequency modes because of quantum effects (fluctuation), 
whereas it is difficult to calculate Cqm(0- 

Rather than using the population relaxation rate 1 /Ti , we could compute the rate of tran- 
sition between pairs of vibrational quantum states 

2n 



where n is the vibrational quantum number. In the limit that ptioJs 3> 1 as we consider 
here, the splitting between vibrational levels is large compared with the thermal energy. At 
equilibrium, the system oscillator will be ground state dominated, and we find that 

1 2Cqm(^5) _ , qm 



kfZo- (4) 



For such a system, we are free to consider the rate of vibrational relaxation in terms of the 
ensemble averaged relaxation rate 1/Ti or the microscopic rate constant fc^^g — results 
will be equivalent. 

In the limit that PhuJs — > 0, on the other hand, the splitting between states becomes much 
smaller than the thermal energy and the results are not equivalent. The rate constant kfT^Q 
diverges, while the population relaxation rate 1/Ti is well behaved 

^ - CqmiuJs)- (5) 

In this work, we will present our results in terms of 1/Ti. 
3.2 Quantum correction factor 

While Cqm(i) is difficult to compute for all but the simplest systems, it is often possible to 
compute the classical analog 

Cci(t) = — (-^(WO))ci (6) 

ms 

for highly non-linear systems consisting of thousands of atoms. The above bracket denotes 
a classical ensemble average. The challenge is to relate the quantum mechanical correlation 
function to its classical analog. An approach explored by Skinner and coworkers has proved 
to be quite productive ^2]- It involves relating the spectral density of the quantum system to 
that of the analogous classical system as 

Cqm(c^s) = Q{(^s)Ccl{^s) (7) 
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detailed balance Q{uj) = Q{—io)e^^'^ and satisfy the "classical" limit that as (3huj becomes 
small, the QCF approaches unity. Using this result, we may rewrite Eq. as 

' -"^^"'kM^s) (8) 



Note that the classical VER rate is defined as = Cdi^s)- 

The QCF for a one phonon relaxation mechanism is 

Qnio^) = i_,-,n. - (9) 

However, as the CD streching mode falls in the transparent region of the DOS (Fig. [2j, a 1:1 
Fermi resonance (linear resonance) is not the possible mechanism of VER. As such, the lowest 
order mechanism available for the VER of the CD mode should involve two phonons. 

We have employed Skinner's QCF approach for two-phonon relaxation If the assumed 
two-phonon mechanism assumes that two lower frequency bath modes, having frequencies uja 
and LOs — UJA, are each excited by one quantum of vibrational energy, the appropriate QCF is 

Qhh{^^s) = Qh{^^a)Qh{ujs - coa)- (10) 

Alternatively, if the assumed two phonon mechanism is one that leads to the excitation of one 
bath vibrational mode of frequency coa, with the remaining energy h{ujs — oja) being transferred 
to lower frequency bath rotational and translational modes, the appropriate QCF is 



Qh-hs{u^s) = QH(u;A)VQ^K^^^^^^e^''^"^""^^/', (11) 

The functions Qh, Qhh, Qh-hs are called the harmonic, harmonic-harmonic, harmonic-harmonic- 
Schofield QCF, respectively. 



3.3 Normal mode calculations for cytochrome c 

To compute the QCF requires a knowledge of, or guess at, the mechanism of vibrational energy 
relaxation. Likely bath modes must be identified and a combination of those modes must meet 
a resonance condition enforced by the conservation of energy in the transition. The candidate 
modes are identified using quenched normal mode (QNM) or instantaneous normal mode (INM) 
calculations. 

In Fig. [21 we show the density of states (DOS) for cyt c in vacuum and in water at 300K. 
These are the INM spectra, so they contain some negative (actually imaginary) components. 
The basic structure of this DOS is similar to that of other proteins like myoglobin |lUlll6j . The 
librational and torsional motions are embedded in lower frequency regions below 2000 cm~^, 
and vibrational motions are located in higher frequency regions around 3000 cm~^. There is 
a transparent region between 2000 cm~^ and 3000 cm~^; the peak due to the CD mode falls 
in this region near 2200 cm~^. The VER of this CD mode is our target in this study. Note, 
furthermore, that the spectra in vacuum and in water are very similar: This indicates that 
water solvent might not affect the simulation results. This conjecture will be confirmed later. 



3.4 Application to VER of the CD bond in cytochrome c 

Bu and Straub employed the QCF approach to estimate the VER rate of a CD bond in 
the terminal methyl group of MetSO in cyt c (Fig. ^ . Their calculations were done using the 
program CHARMM and cyt c was surrounded by water molecules at 300K. Compared to 
this, we have used molecular dynamics simulations of cyt c in vacuum at 300K to compute the 
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Figure 2: Density of states p{uj) for cytochrome c in vacuum (solid line) and in water (dashed 
line) at 300K calculated by INM analysis. 

classical autocorrelation function for the force acting on the same CD bond. The results have 
been used to make estimates of both l/Tf^ and 1/T^'~'^ . 

In Fig. El the force autocorrelation function and its power spectrum are shown for four 
different trajectories. We have observed that the force fluctuation and the magnitude of the 
power spectrum for cyt c in vacuum is very similar to those computed for cyt c in water. We 
conclude that the effects of water on the VER rate are negligible. With the CD bond frequency 
u>s = 2133 cm~^, we find 1/Tf = Cci('i^5) — 1 ps~^, i.e. the classical VER time is about 1 ps. 

To apply QCFs for two-phonon relaxation, Eqs. (jTU)) and ((TT|) . to this situation, we need 
to know UJA- We have found that the CD mode is strongly resonant with two lower frequency 
modes, 1655th (685.48 cm^-*^) and 3823rd (1443.54 cm"-*^) modes because \ujs — ^1655 —'^38231 = 
0.03 cm~^. We might be able to choose loa = 1443.54 cm~^ or 685.48 cm~^. 

In Fig. |11 we show uja dependence of the normalized QCF, i.e. Q = Q/ {/Shoos) = Tf'/T^^^^ 
at 300K and at 15K. If we choose loa = 1443.54 cm~^ at 300K, Q = 2.3 for the harmonic- 
harmonic QCF and 2.8 for the harmonic-harmonic-Schofield QCF. Thus we have T^*^^ = 
Tf/Q = 0.3 ~ 0.4 ps. It is interesting to note Q at 15K varies significantly depending on the 
QCF employed. We will discuss this feature later. 

3.5 Fluctuation of the CD bond frequency 

We have discussed the fluctuation of the frequency for the CD bond |11,. In the equilibrium 
simulation, the instantaneous normal mode analysis has been employed for each instant of time 
to generate a time series iocoit) for the CD bond frequency. From this time series, we can 
calculate the frequency autocorrelation function 

1 fT 

C{t) = 5oocD {t)5LVcD (0) = 7^ / drSiOcD {t + t)6ujcd (r) (12) 

J JO 
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Figure 3: Left: Classical data for the force-force correlation function. Middle: Fourier spectra 

for the correlation function. Right: Magnification of the middle figures around the CD bond 
frequency. These data are taken from four different trajectories of the equilibrium simulation. 
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Figure 4: Normalized harmonic-harmonic (HH) and harmonic-harmonic-Schofield (H-HS) QCF 
at 300K (left) and at 15K (right). 



where the overline means a long time (T) average, and dujcnit) = uJcoit) — uJcDif)- The 
correlation time is defined as 

where (Aw)^ = C(0). From Fig. El we found Aa; ~ 8.5 cm^^ and Tc — 0.2 ps. Since AwTc <C 1, 
according to Kubo's analysis Qf2n-, the lineshape should be homogeneously broadened, i.e., its 
shape is Lorentzian. This is also the case for cyt c in water We also confirmed that the 
potential barrier of the methyl group to rotate is significantly greater than the thermal energy 
(barrier height ~ 3 kcal/mol > thermal energy ~ 0.6 kcal/mol) so that we do not expect 
inhomogeneity in the line shape. 

These results support the validity of employing a normal mode type study of VER in cyt c 
as the structure of cyt c is rather rigid around the CD bond and the dynamics of the bond, on 
the scale of VER, should be well modeled by NMs. Of course, to describe VER among NMs, 
we must include nonlinear coupling terms. In the next section, we will discuss this reduced 
model approach for VER in cyt c. 



4 Reduced model approach 

The QCF approach is attractive in that it allows us to base our estimate of the VER rate 
on a force autocorrelation function that contains the non-linearity of the coupling between 
the system and bath modes to all orders. To obtain an estimate of the VER based on a 
more accurate representation of the system's quantum dynamics, we expand the potential as 
a Taylor series in the normal modes of the of the system and bath. In this representation, 
the system and bath coordinates are "normal modes;" to second order in the expansion of 
the potential energy, the system and bath modes are uncoupled and non-interacting. The 
interaction "coupling" between the system and bath modes first appears at third-order. In this 
section, we describe a perturbation theory estimate of the rate of VER of the CD mode that 
represents the system-bath coupling to lowest, third-order in the system and bath coordinates. 
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Figure 5: The frequency autocorrelation function. The data up to 20 ps were used to calculate 
the correlation function. 



4.1 Reduced model for a protein 

The reduced model approach utilizes the normal mode picture of a protein, expanding the 
residual term perturbatively as 
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Thus the force applied to the system mode is 



^ = = ~3 ^ Gs,k,iqkqi - 4 X! Hs,k,i,mqkqiqm + 

k,l k,Lm 



(19) 



where we have used the permutation symmetry of Gkim and Hkimn- If it is enough to include 
the lowest order terms proportional to Gkirm substituting them into Fermi's golden rule Eq. 
we can derive an approximate VER rate as ^HI 
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following 



Cff = —-^^{l + nk + ni + 2nkni), (21) 

Ci7^ = ^^^K + n; + 2nfcn0, (22) 

^kl = -^Gs,k,u (23) 

= l/{e^^^^-l). (24) 

4.2 Maradudin-Fein formula 

There exists another well known formula to describe the VER rate, the Maradudin-Fein (MF) 
formula P|ITlj. 



W = W^decay + W^colb (25) 

n 



W^coii = 7 ^ (27) 

with a width parameter 7. Note that Eq. ()2U() and Eq. (|25() are equivalent in the limit of 7 — > 
as shown by Kenkre, Tokmakoff, and Fayer p^]- However, they disagree with a finite width 
parameter such as 7 ~ 100 cm~^. In this chapter, we use the MF formula and consider its 
classical limit {Ti 0) defined as 

{A^i? ( 1 , 1 



wcl ^ ^ V ^ ^''^ _ _L ± i C281 

"^^''^ 2msPoosf^ u^k^l W ^iJ I' + iu^S-u^k-u^l)^' ^ ' 

We note some properties of the formula: Wdecay > l^dLay ^^'^ W^coll < ^coll which is derived 
from 

l/(e^ - 1) - l/{ey - 1) 
1/x - 1/y 

for X, y > 0. We can define an effective QCF as 



< 1 (31) 



This should be compared with the normalized QCFs [Q = Q{uJs)/{f3huJs)] found in the litera- 
ture. 
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To apply the Maradudin-Fein theory to the case of CD bond relaxation in cytochrome c, we 
must numerically compute the third order coupling elements. As the protein has more than 
10^ modes, there are more than 10^ third order coupling elements representing the coupling of 
the "system" CD bond to the "bath" modes of the surrounding protein and solvent. 
We have employed the finite difference approximation 

2dqsdqA~ 2^^^"^^^" 2Aqs ^^^^ 

where Uik is an orthogonal matrix that diagonalizes the (mass-weighted) hessian matrix at the 
mechanically stable structure Kij, and Kij{zizAqs) is a hessian matrix calculated at a shifted 
structure along the direction of a selected mode with a shift ±Aqs. 

Note that, in the large number of coupling elements, most are small in magnitude. Of those 
that are larger, most fail to meet the resonance condition and do not contribute significantly 
to the perturbative estimate of the VER rate. See for the details. 



4.4 Width parameter 

We show the width parameter 7 dependence of the VER rate in Fig. |HJ^ We consider other 
lower frequency modes (wsaso = 1330.9 cm"^, wigge = 829.9 cm^^, ujie^^ = 685.5 cm~^) as 
well as the CD mode {cocD = 2129.1 cm^^) for comparison. From the former analysis of the 
frequency autocorrelation function Eq. ()12|). we might be able to take 7 ~ Aa; ~ 3 cm~^ for 
the CD mode, and we have Ti ~ 0.2 ps, which agrees with the previous result with QCFs: 
T^^^ = 0.3 ~ 0.4 ps. 

We also see that the lower frequency modes have longer VER time, a few ps, which agrees 
with the calculations by Leitner's group employing the MF formula ^Sj. The main contribution 
to the VER rate at 7 = 3 cm~^ comes from 1655th mode (685.5 cm~^), a heme torsion, and the 
3823rd (1443.5 cm~^) mode, an angle bend in Met80 (~ 20%). Interestingly, we can conceive a 
peak around 7 = 0.03 cm~^. Given this width parameter, the contribution from the two modes 
is more than 90%. In any case, we can say that 1655th and 3823rd modes are resonant with 
the CD mode because they satisfy the resonant condition (|a;i655 +^3823 — ^cd\ — 0.03 cm~^) 
and the coupling elements between them is relatively large (1^16553823! — ^-^ kcal/mol/A'^). 

However, this close resonance does not necessarily lead to the conclusion that it forms the 
dominant channel for VER of the CD strech. There is a competing near-resonance involving 
the 3330th mode (1330.9 cm~^), an angle bending mode of Met80, and the 1996th mode (829.9 
cm~^), a stretch-bend mode in Met80. While the resonance is not close (it is within 31.7 cm~^), 
the coupling element is quite large (1^3330 iggel — ^2. 3 kcal/mol/A^). With a larger value of 
7 = 30cm~^, this combination of bath modes becomes the dominant channel for VER of the CD 
stretch. Clearly, the uncertainty in our force field, used to compute the vibrational frequencies, 
and the value of 7, which is rather poorly defined, prevents us from concluding that one or 
another of these two channels will dominate VER of the CD stretch at room temperature. 

In the left of Fig. d we show the effective QCF calculated from Eq. (|32|) at 300K, which 
is Qeff — 2.3 for the CD mode with 7 = 3 cm~^. This value better agrees with the (normal- 
ized) harmonic-harmonic QCF [Eq. (|10j) ]. compared to the harmonic-harmonic-Schofield QCF 
[Eq. pij) ]. The Qeff for the other modes are more or less unity, which indicates that these 
modes behave classically at 300K.'^ In contrast, as is shown in the right of Fig. [3 Qeflf at 15K 

^Note two limiting cases of 7 dependence: l/Ti oc 7 when 7 is very small, and l/Ti oc I/7 when 7 is very 
large. This is easily recognized from the Lorentzian form. 

^We notice an interesting behavior for 1655th mode, i.e. Qeff becomes very much smaller than unity at 
7 ~ 0.03 cm"^. In this case, we observe that Wcoii ^ VKdecay because of the resonance: ojiess +ci;3823 —ujcd — 
(actually 0.03 cm~^). In such a case, Qcs becomes less than unity because VKcoii < W^con- 
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Figure 6: VER rates for the CD mode {ujcd = 2129.1 cm ^) and the other lower frequency 
modes (waaao = 1330.9 cm^^, wigge = 829.9 cm~^, wiess = 685.5 cm"-*^) as a function of 7 at 
300K (left) and at 15K (right). 
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Figure 7: Effective QCF for the CD mode {ujcd = 2129.1 cm ^) and the other lower frequency 
modes ((^3330 = 1330.9 cm~^, cjiggg = 829.9 cm~^, wiess = 685.5 cm"-^) as a function of 7 at 
300K (left), and at 15K (right). 
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it is proportional to the temperature (see also Fig. IH)). A similar trend is found in the right 
of Fig. ^ where the harmonic-harmonic QCF (Q ~ 40) is comparable to QeS- On the other 
hand, the harmonic-harmonic-Schofield QCF gives an exponentially large value of Q, showing 
strong deviations from Qeflf- We should bear in mind that different QCFs lead to significantly 
different conclusions at low temperatures. 

4.5 Temperature dependence 

In Fig. IHl we show the temperature dependence of the quantum and classical VER rate cal- 
culated by the MF formula and the classical limit of the MF formula. At high tempratures 
(~ lOOOK), the quantum VER rate agrees with the classical one, but they deviate at low tem- 
peratures. The former becomes constant due to the remaining quantum fluctuation (zero point 
energy) whereas the latter decreases as oc (temperature). The "cross over temperature" where 
the VER behaves classically is smaller for the lower frequency modes compared to that of the 
CD mode as expected. 




Figure 8: Quantum (left) and classical (right) VER rates for the CD mode {ujcd = 2129.1 
cm~^) and the other lower frequency modes (wasao = 1330.9 cm^^, wiggg = 829.9 cm~^, 
"^1655 = 685.5 cm~^) as a function of temperature with 7 = 3 cm^^. 



5 Discussion 

5.1 Comparison with experiment 

Here we compare our results with the experiment by Romeberg's group JH]- They measured 
the shifts and widths of the spectra for different forms of cyt c; the widths of the spectra 
(FWHM) were found to be AwpwHM — 6.0 ~ 13.0 cm^^. From the discussions of Sec. 13. 5| we 
can theoretically neglect inhomogeneous effects, and estimate the VER rate simply as 

Ti ~ 5.3/Aa;FWHM (ps) (34) 

which corresponds to Ti ~ 0.4 ~ 0.9 ps. This estimate is similar to the QCF prediction 
using Eq. Q (0.3 ~ 0.4 ps) and the reduced model approach using Eq. p0|) or p5|l (0.2 ~ 
0.3 ps). We note that this value should be compared to the VER time of the C-H stretch 
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studies, e.g. on temperature dependence of absorption spectra or on time-resolved spectroscopy, 
will clarify which methodology is more applicable. 

Romesberg's group studied Met80-3D (methionine with three deuteriums) while we have 
examined Met80-1D (methionine with one deuterium). In the case of Met80-3D, there are three 
peaks in the transparent region. It should be possible to consider the VER of each of the three 
modes, to make a more direct comparison of the predictions of our theoretical models with the 
results of their experimental studies. 

5.2 Validity of Fermi's golden rule 

We next discuss the validity of our approaches. Since our starting point is the perturbative 
Fermi's golden rule, our two approaches should have a limited range of validity. Naively speak- 
ing, the force applied on the CD mode should be small enough, but how small it should be? 

We follow Kubo's derivation of a quantum master equation using the projection operator 
technique j22j. He derived an equation for the evolution of the system density operator a(t) 

^fii) = f_^dT[q{t)q{T)a{T)^{t-T)-q{t)a{T)q{T)^{-t + T) 

+a{T)q{T)q{t)^{-t + r) - q{T)a{T)q{t)^{t - r)]. (35) 

The interaction Hamiltonian is assumed to be T^int = as in our case, i.e. q is the 

system coordinate and J- mainly contains the bath coordinates. We have defined the force 
autocorrelation function $(t) = TrB{pBJ^(t)J^(0)} = (j^(t)J^(O)). Note that Eq. (jSSJ) is just a 
von Neumann equation using the projection operator technique, and it is not a master equation 
yet. 

If ^{t) decays fast, we can replace <j{t) in the integral with cr(t), and the dynamics be- 
comes an approximate Markovian dynamics. If this approximation is valid, Fermi's golden rule 
describes the relaxation dynamics of a{t) |22j. The validity of the golden rule relies on the 
validity of the Markov approximation. 

From Eq. ()35|) . the relaxation rate of (T{t) can be estimated as 

1/r, ~ ((g2)^V^')Tc (36) 

where we have assumed 

$(t) ~ .F2e-l*l/^^ (37) 
The Markov approximation [o"(r) ~ 0"(i)] holds for 

> Tc. (38) 

We have a criterion for the validity of the Markov approximation 

e = {i)T'^TllT? <. 1. (39) 

In our case as well as the case of HOD in D2O, the ratio is "just" small (see Table Applying 
Fermi's golden rule to these situations should be regarded as a reasonable estimate of the VER 
rate. As alternative approaches that avoid this underlying Markov approximation, one can 
employ more "advanced" methods as mentioned in Sec. IHl 

5.3 Higher order coupling terms 

Up to now, we have only included the third order coupling terms to describe the VER of 
the CD mode. However, we must be concerned with the relative contribution of higher order 
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Table 1: The parameters in Eq. H39|) for various molecules in AKMA units (unit length = 1 A, 
unit time = 0.04888 ps, unit energy = 1 kcal/mol). The data for HOD in D2O, CN~ in water, 
and CO in Mb are taken from [^Tj, and respectively. 











e 


CD in cyt c 


0.01 


5.0 


1.0 


0.5 


HOD in D2O 


0.01 


5.0 


1.0 


0.5 


CN~ in water 


0.002 


1.0 


0.5 


0.01 


CO in Mb 


0.002 


1.0 


1.0 


0.02 



mechanism, e.g. the contribution due to the fourth order coupling terms in Eq. (|18|) . This 
is a very difficult question. As there are many terms (~ 10^) included, we cannot directly 
calculate all of them for cyt c. We have found that it is not sufficient to include only the third 
order coupling terms to reproduce the fluctuation of the force on the CD bond. However, this 
does not necessarily mean that the VER rate calculated from the third order coupling terms is 
inadequate. 

The main contribution from the fourth order coupling terms to the VER rate and the force 
fluctuation are written as 

^ / 1 \ ^ \Hs,k,l,m\' ^^^^ ^^)^ (40) 



and 



A{5J^^) ~ , (41) 



|2 



respectively. Even if A{6J-'^) becomes large, A(l/Ti) is not necessarily large because of the 
resonance condition {ujs — ujk — ijJi — ojm ^0). It is a future task how to evaluate the effects 
due to higher order coupling terms. It might be interesting to compare the classical limit 
using the reduced model approach and the Landau- Teller-Zwanzig approach because there is 
no ambiguity about how to choose QCFs 



6 Summary 

In this chapter, we have examined VER in a protein from the QCF approach and the reduced 
model approach, and compared the results. For the CD mode in cyt c (in vacuum) at room 
temperature, both approaches yield the same result for the VER rate, which is also very similar 
to an estimate based on an experiment by Romesberg's group. Our work demonstrates both 
the feasibility and accuracy of a number of theoretical approaches to estimate VER rates of 
selected modes in proteins. 

The QCF approach is appealing in that the calculation of the force autocorrelation function 
is straightforward and feasible, even for systems of thousands of degrees of freedom. Moreover, 
the classical force autocorrelation function includes all orders of non-linearity in the interaction 
between the system oscillator and the surrounding bath. A weakness of the QCF approach 
is that we don't know which QCF to choose a priori. We must assume a mechanism for 
VER before computing the rate. Moreover, the temperature dependence of the rate of VER 
is sensitive to the mechanism, whether it involves few phonons or many phonons. The choice 
of the form of the quantum correction factor can make a significant difference in the predicted 
rate of VER at lower temperatures. 
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of the reduced system is accurately treated. Using the reduced model approach, there are 
two ways to estimate the quantum mechanical force autocorrelation function: (1) numerical 
calculation of the quantum dynamics for a model Hamiltonian of a few degrees of freedom, 
including all orders of non-linearity in the potential, and (2) analytical solution for the quantum 
dynamics using perturbation theory that includes many bath modes but only the lowest order 
non-linear coupling between system and bath modes. We have employed the latter approach 
through the use of the Maradudin-Fein formula. A weakness of our reduced model approach is 
that the method neglects the higher order coupling elements beyond third order, which cannot 
be justified a priori (301 1^- 

We have pursued a comparative study in which we seek consensus in the estimates of 1/Ti 
that result of the QCF approach and the perturbation theory. A rather remarkable result of our 
study is that while the absolute value of the quantum corrections to the classical VER theory 
are large (on the order of a factor of 40), the results of the QCF and the perturbation theory 
approaches are in close agreement. This is all the more remarkable given the fact that the results 
of the perturbation theory require a calculation of the third-order coupling constants and the 
estimation of the "lifetime" parameter 7. As we have shown, the dominant channel for VER, 
derived from the perturbation theory, depends upon the choice of 7. For smaller 7 = 3cm~^, 
the dominant mechanism is a close resonance (within 0.1 cm~^), a combination of a heme 
torsion and MetSO angle bending mode, with a weak coupling. For larger 7 = 30cm~^, the 
dominant mechanism appears to be a less perfect resonance (within 31.7 cm~^), a combination 
of a different angle bending mode and a bend-stretch mode in Met 80, with a strong coupling. 
Such detailed knowledge of 7 is essential to predict a mechanism for VER. 

Our study raises two important questions: (1) What is the optimal set of coordinates for 
modeling and interpreting VER in proteins? In the QCF approach, we treated the relaxing 
bond (CD bond) as a local mode that is coupled to vibrational modes of the bath. In the 
reduced model approach, on the other hand, we treated all the vibrational modes including the 
relaxing mode (CD mode) as normal modes that are coupled each other with the third order 
nonlinear coupling terms. Our numerical results showed that the two approaches give similar 
results for the VER rate of the CD bond or mode, but it remains to be seen that this is a kind 
of coincidence or there is a theoretical ground of their equivalence (if the QCF is appropriately 
chosen). (2) What is the physical origin of the width parameter 7 and how to calculate it? 
In this work, we suggested to use the relation 7 ~ Auj where Aa; represents the fluctuation 
of the CD mode (or bond) frequency. We think this is reasonable but there is no theoretical 
explanation of this. If the VER rate does not significantly depend on 7, this is not a serious 
problem, but this is not always the case. Thus we need an "ab initio" way to derive the width 
parameter 7. One appealing way is to regard 7 as a hopping rate between potential basins 
(inherent structures) 32, 33_j- The other is the more rigorous quantum mechanical treatment 
of the tier structure of energy levels in the protein j34j . 

The results of our study are derived through the use of an approximate empirical energy 
function (force field) which has not been "tuned" to provide accurate frequencies of vibration for 
all protein modes. Our predicted rates of VER depend sensitively on the closeness of the reso- 
nance between the system and bath modes. Clearly, we must resort to the re-parameterization 
of the empirical potential to fit with experimental data or higher levels of theory (ab initio 
quantum chemistry calculation) in an effort to refine our estimates of the frequencies of vibra- 
tion and the details of non-linear coupling between vibrational modes of the protein. This is a 
challenge for both experimental and theoretical studies. 

Recent advances in experiment and theory make the present time an exciting one for the 
detailed study of protein dynamics. A variety of methods have been applied to examine VER 
in molecules, including nonequilibrium MD methods [321, time-dependent self-consistend field 
methods |311 136j . mixed quantum-classical methods ^T, and semiclassical methods [381 139| . 
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or 2D-IR signals j4Ul I41j as probes of protein structure and dynamics. Extensions of these 
studies will provide us with an increasingly detailed picture of the dynamics of proteins, and 
its relation to structure and function. 
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